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Abstract 

We study algebraic algorithms for expressing the number of non-negative 
integer solutions to a unimodular system of linear equations as a function 
of the right hand side. Our methods include Todd classes of toric varieties 
via Grobner bases, and rational generating functions as in Barvinok's 
algorithm. We report polyhedral and computational results for two special 
cases: counting contingency tables and Kostant's partition function. 



1 Introduction 

The object of study in this paper is the vector partition function 

— # { a; : Ax — b,x > 0, x integral }, 

where ^ is a fixed d x n-unimodular integer matrix and 6 is a variable vector in 
Z'K Here we say that A is unimodular if the polyhedron {x : Ax = b, x > 0} has 
only integral vertices whenever b is in the lattice spanned by the columns of A. 
This is a slight generalization of the definition of "unimodular" used in §19]. 
We further assume that Ker{A) n M"g — 0, which is equivalent to 0a (^) < oo 
for all b. We regard (/)a as a function on cone(^), the cone of non- negative linear 
combinations of the columns of A, since 0a (^) = if 6 is not in cone (A). The 
following result about vector partition functions is well-known (see e.g. |23|| ): 

Theorem 1.1 The function 4>a is piecewise polynomial of degree n — rank(A). 
Its domains of polynomiality are convex polyhedral cones, called chambevs of A. 



The main purpose of this paper is to develop practical methods for uni- 
modular counting. By unimodular counting we mean preprocessing the given 
unimodular matrix A and generating the polynomials for 0a on the various 
chambers. Each output polynomial is represented either explicitly as a sum of 
monomials, or implicitly as an oracle which allows for quick evaluation of 0a 
at any b in that chamber. Unimodular counting has many applications, rang- 
ing from statistics and randomized algorithms pq| to representation theory 
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[ p^l PH . For instance, the widely known problem of counting contingency tables 
is the case when A is the incidence matrix of a complete bipartite graph [ill . 

Our benchmark on unimodular counting is the work of Mount . His 

approach is based on interpolating the chamber polynomials, by evaluating 
4>A{b) for sufficiently many right hand sides 6, coupled with divide- and-conquer 
decompositions and advanced parallel computation techniques. Both the eval- 
uation and the divide-and-conqucr schemes depend on the specific matrix A. 
Mount reports the complete solution for contingency tables of size 4x4. In 
Welsh's survey on approximate counting. Mount's computations for 4x4- 
tables are mentioned as the limit for exact counting on today's computers. 

Mount's method does not take full advantage of the rich algebraic structure 
underlying (j)A ■ On page 64 of his thesis ||l8| , he writes "There are some results in 
commutative algebra that relate the (chamber) polynomials to "Hilbert series" 
and "Todd classes", but these structures encode a lot of information and are in 
themselves hard to compute. The strategy taken here is to assume access to a 
counting oracle .... and then recover the desired polynomial by interpolating..." 

We shall demonstrate that algebraic algorithms perform much better than 
Mount had surmised. In fact, using rather simple test implementations, we can 
now count 4x5 and 5x5 contingency tables with arbitrarily large margins. 

The algebraic methods apply to any unimodular matrix A and work in- 
dependently of the size of the right-hand-side vector b. In fact, our original 
motivation for this project was the open problem, stated by Kirillov page 
57], of computing the number of chambers for Kostant's partition function of 
the root system v4„i_i. In this special case, our unimodular matrix is the inci- 
dence matrix of the complete directed graph K„i. We solve Kirillov's problem 
for m < 7, and we compute all chamber polynomials up to to < 6. Using these 
polynomials we provide an on-line calculator for Kostant's partition function at 



www.math.ucdavis.edu/~deloera/kostcLnt.html. We also prove some other 



new results on the geometry of chamber complexes of unimodular matrices. 
This paper presents two algebraic algorithms for unimodular counting: 

1. A Grobner bases algorithm, which computes the Todd class of the toric 
variety defined by our polytope {x >Q : Ax = 6}, is given in Section 2. 

2. The BBKLP method, which computes the generating function for all lat- 
tice points in the polytope {x > : Ax — b}, is given in Section 3. 

The acronym BBKLP refers to five mathematicians: Barvinok, Brion, Kho- 
vanskii, Lawrence, and Pukhlikov. The most important complexity result in 
our area is Barvinok's polynomial-time algorithm for counting lattice points in 
rational polytopes of fixed dimension |9|. Barvinok's algorithm is based on 
earlier work by Brion, Khovanskii, Lawrence, and Pukhlikov. For a complete 
bibliography see the survey article of Barvinok and Pommcrsheim [Q. When 
A is unimodular, Barvinok's algorithm specializes to the BBKLP method and 
runs very fast in practice. This answers a question of Mount p8| , page 56]. 
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We implemented methods (1) and (2) in the computer algebra packages 
Macaulay 2 and Maple respectively. Details are described in Sections 4 and 5. 
We expect a significant further speed-up by combining our algebraic approach 
with Mount's parallel computing techniques. In a future project we will extend 
the various methods for computing (/>^ to non-unimodular matrices A. 



2 Method One: Counting Using Grobner bases 

We describe now our first algebraic algorithm for solving the following counting 
problem associated with any unimodular dx n-matrix A: Determine the number 
4>A{b) of non-negative integer solutions u £ N" of the hnear equations A-u = b. 

The following discussion makes use of well-known facts from algebraic ge- 
ometry (see 0); specifically, we demonstrate how to effectively compute the 
Todd cohomology class of a toric manifold defined by a unimodular matrix. 

Our running example is the following unimodular 3 x 5-matrix: 

/l 1 l\ 
A = 10 10. 

\0 1 1/ 

The vector partition function for this matrix equals 

{bc + b + c + 1 if a > b + c and b,c> 0, 

+ fa + 1 if min{6, c} > a > 0, 

ab-\b^ + \b + a + l if c > a > 6 > 0, 

ac — ^c? + ic + a + l ii b > a > c, 

a6+ac-i(a^+f)^ + c^) + i(a+6+c) + l if 6 + c > a > max{6,c}. 

For our exposition it is more convenient to express the vector partition function 
as i>A ■ N" — *■ N where iPa{v) is the number of solutions u e N" to the equation 
Au = Av. Clearly, ipA and (pA are related by a simple transformation. For 
instance, in our example we have 'ipA{a, b,c,d,e) = (j)A{a + d -\- e,b + d,c + e). 

The chamber complex of a unimodular matrix A is defined as the common 
refinement of all triangulations of A. For the 3 x 5-matrix A above, the chamber 
complex is the given subdivision of cone(A) = R>g into five triangular cones. 
We refer to and for details on chamber complexes and ||2^ for an intro- 
duction to Gale transforms and triangulations. We assume that rank(j4) = d. 

Lemma 2.1 The chambers of A are in bijection with the regular triangulations 
of any Gale transform A of A. Non-regular triangulations of A are in bijection 
with the virtual chambers of A. 

Thus generating the chambers of our unimodular matrix A is the same as 
generating all regular triangulations of a Gale transform A. It is well-known 
that the regular triangulations can be generated by a applying bistellar flips to 
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a seed regular triangulation (see |^). Bistellar flips are topological operations 
that transform a triangulation into another. One has to be careful as sometimes 
a flip creates non-regular triangulations, but regularity of a triangulations can 
be checked by linear programming. When necessary we have performed these 
calculations using the software packages Puntos Q and Topcom pO| . 

We first characterize the chamber complex in algebraic terms. Let S = 
k[xi, . . . , Xn] be the polynomial ring over a field k which contains the rational 
numbers. The variables of S index the columns of the matrix A — (aij). Let 
Ja denote the ideal in S generated by the binomials x^^Xj'^ • ■ • a;"*" — 1 for 
i = 1,2, ... ,d. For any positive weight vector w G R", let in^^{JA) denote the 
ideal generated by the ui-initial forms of the binomials in J a. li w is generic, 
then in^{JA) is a monomial ideal. It was shown in Corollary 8.9] that the 
matrix A is unimodular if and only if all initial monomial ideals in^^JA) are 
square-free. Two weight vectors w and w' in K" lie in the same cone of the 
Grobner fan if inw{JA) — in^'{JA)- By the results in p4| §8] this happens if 
and only if, for every linearly independent subset a — {ai-^, . . . , a^^} of column 
vectors of A, the vector Aw lies in the cone spanned by a if and only if the 
vector Aw' lies in the cone spanned by a. This implies the following result: 

Proposition 2.2 The chamber complex of A equals the Grobner fan of Ja- 

Algebraic algorithms for computing Grobner fans are described in |2^, §3]. 
The state of the art on this subject is the work of Huber and Thomas pGf . 
We now explain how to compute the polynomial representing ^a on any given 
chamber. Suppose that w is a positive integer vector in the interior of that 
chamber. Then M = inw{JA) is a square- free monomial ideal. It was shown 
in Corollary 7.4] that M encodes the face poset of the simple polytope 

= {u : u> and Au^ Aw}. 

For any (n— d)-element subset / of {1, ... , n}, the equations Ui = 0, i € I define 
a vertex of P^, if and only if {xj : j ^ /) is a minimal prime of M . Writing S^j 
for the normal fan of the simple polytope P^,, this can be restated as follows: 

Proposition 2.3 The Stanley- Reisner ideal of the fan equals M — inm{JA). 

In our running example, with w — (1,1,1.1,1), the polytope Pyj is a pen- 
tagon and the fan has five rays in the plane. This is encoded by the ideal 

M ^ {A,B,C) n {A,B,E) n {B,D,E) n {C,D,E) n {A,C,D). (1) 

Returning to the general case, our goal is to count the lattice points in the 
polytope Pw We use known methods from toric geometry for this computation. 
An introduction can be found in Section 5.3 in Fulton's book See also Q. 

Let denote the projective toric variety defined by the fan E^. The variety 
Xw is smooth, for all generic w, since A is unimodular. Let La denote the ideal 
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in S ~ k[xi, . . . , x„] generated by the linear forms bixi + • • • + &ri.T„, where 
(bi, . . . ,bn) runs over all vectors in the kernel of the matrix A. The cohomology 
ring of Xyj with coefHcients in our field k is the artinian graded fc-algebra 



Arithmetic operations in this algebra are performed using normal form reduction 
relative to any Grobner basis of the ideal M + L^. Since is an irreducible 
complex manifold of dimension n — d, the top cohomology group H^"~^'^{X^, k) 
is a one-dimensional vector space. There is a canonical choice of a basis vector for 
that one-dimensional fc-vector space, namely any square-free monomial Wi^j Xi 
which indexes a vertex of Pyj. This is equivalent to {xj j ^ I) being a minimal 
prime of M. Since X^ is smooth, any two such monomials are congruent to 
each other modulo M + La- The resulting element of H*{Xyj;k) represents 
the cohomology class which is Poincarc dual to a point. 

The following rule uniquely defines a /c-linear functional called the integral: 



Writing top(p) for the degree n — d component of p, we require that top{p) — 
(/ P) ' Yliei ^^^^ in M + La, where / is any index set as above. 
Algorithm 1. (Computing the integral of a cohomology class of X^) 
Input: A polynomial p{xi, . . . , a;„) with coefficients in a field k D Q 
Output: The integral p of the corresponding cohomology class on X^. 

1. Compute any Grobner basis Q for the ideal M + La- 

2. Let m denote the unique standard monomial of degree n — d- 

3. Find any minimal prime {xj : j ^ I) of M, and compute the normal form 
of riie/ modulo the Grobner basis Q- It looks like 7 • m, where 7 is a 

non-zero element of k- 

4. Compute the normal form of p modulo the Grobner basis Q, and \ei 5 & k 
be the coefficient of m in that normal form. 

5. Output the scalar 5/7 G k. 

To compute the number of lattice points in P^, we note that there is a 
special element in the cohomology ring i?*(X^; fc), denoted td{xi, . . . , a;„) and 
called the Todd class of the toric variety Xyj- The Todd class is represented 
(non- uniquely) by a (non- homogeneous) polynomial with rational coefficients in 
the variables xi, - . - ,Xn- The polynomial td{xi, . . . , a;„) does what we want: 



H*{X^;k) = 0ij2'^(X^,fc) = S/{M + La)- 



(2) 



r=0 
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Here the exponential of a linear form in (|2[) is defined by the terminating series 

n n — d^ 
eXpC^WiXi) = — ■ {wiXl + W2X2 + \-WnXnY- (3) 

Pommersheim [Q gives an algorithm for computing the Todd class, which works 
efficiently even for non-unimodular A. For our applications, however, we prefer 
to use the basic formula given in the first line on page 110 in Fulton's book |p^ : 

In this expansion we list only terms of degree < n — d, so that becomes a 
polynomial in xi, . . . , a;„ with Q-coefficients. We conclude with our main result. 

Theorem 2.4 The following algorithm computes the polynomial that represents 
ipA on a chamber containing a given non-negative vector w G M"; 

1. Determine the linear inequalities defining the given chamber. 

2. Let M be the ideal generated by the leading monomials of the Grobner 
basis for J a with respect to w and compute the ideal representing the 
kernel La of A. Use these two ideals to construct the cohomology ring. 

3. Apply Algorithm 1 to the product of the polynomials in (^) and (|^). 

A main advantage of this algorithm over other methods is that we can do 
the computation parametrically, over the field k = Q(it;i, . . . , Our output 
is the actual polynomial for ijjA, not just some numerical evaluation of it. 

For our running example we take the polynomial ring S — k[A^ B, C, D, E] 
over the field k = Q(a, b, c, d, e). We fix the reverse lexicographic Grobner basis 
for the ideal M + La, where La ^ {A + B- D,A + C~E) and M = in^{JA) 
is the monomial ideal in (|l|). The Todd class is computed from the formula 

{l+A/2+A^/l2){l+B/2+B^/12){l+C/2+C^/12){l+D/2+D^/12){l+El2+E^/l2) 

The normal form of this expression with respect to our Grobner basis equals 

td{A,B,C,D,E) = DE + C/2 + D + E/2 + 1. (5) 

Likewise, the exponential of the general divisor (|^) on our toric surface, 

l + {aA + bB + cC + dD + eE) + ^{aA + bB + cC + dD + eEf, 

has the following normal form with respect to our Grobner basis: 

1 + {a~b+e)E+{b+c-a)C+{b+d)D + {ab+be+ac+cd+de-{a^ -b^ -c^)/2)DE 

Multiply this expression with (^, reduce it to normal form, and extract the 
coefficient of the standard monomial DE. The result is the desired polynomial 
that represents ipAio-, b, c, d, e) on the fifth chamber. Now set d = e = 0. 



6 



3 Method Two: BBKLP generating functions 



In the BBKLP method we associate with any rational polyhedron P in R" the 
following rational generating function in n variables: 



fiP,x) = ^ a;" 

where x" denotes . Brion |^ proved the following result: 

Theorem 3.1 For any rational polyhedron P in R", 

/(P,x) = /(cone(P,i;),x), 

vcrticcs(P) 

where cone(P, = {u G M" : v + 5u E P for all sufficiently small (5 > 0}. 

Each of the series /(cone(P, w), a;) is a rational generating function, which 
can be computed using commutative algebra methods (Hilbert series). For us 
the only relevant case is that of unimodular (also called primitive) cones. A 
unimodular cone is a pointed simplicial cone with generators {ui, . . . , Ufe} that 
form a basis for the lattice M{mi, . . . , u^;} n Z". For unimodular cones, the 
rational generating function takes the following simple form: 

/(cone(P,z;),a;) = Ht^- (6) 

If P is a rational convex polytope then /(P, x) is a polynomial, and this polyno- 



mial has a "short" representation as a rational function by Theorem 3.1. Eval 



uating f{P,x) at x = (1, 1, 1, ... , 1) gives the number of integer points in P. 



However, if we are given /(P, x) as a sum of rational functions as in Theorem 3.1 
then this evaluation is a nontrivial problem since the point x = (1,1,1,...,! 
is a pole of (||). We present our solution to this problem in Algorithm 3. 

For a one-dimensional example, let P be the line segment [0, b]. Then 

/(P,x) = + T = l + x + x"^ + ---+x\ 

1 — X 1 — X ^ 

The value of this polynomial at a; = 1 equals b+ I, the number of lattice points 
in the segment, but the substitution x — > 1 must be performed with care. 

Consider now the polytope P^ = {a; € R" : a; > 0, Ax = b}, where A is 
unimodular and b is in the lattice spanned by the columns of A and lies in the 
relative interior of a maximal chamber. Under these hypotheses, Pfc is a simple 
polytope such that cone(Pb, v) is unimodular for every vertex v of Pf,. We shall 
give a combinatorial formula for the rational functions representing these cones. 
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Consider any subset cr C {1, . . . , n} which is a column basis of the matrix 
A, and let Va- denote the unique vector in R" with support a and satisfying 
A ■ Va — h. The entries of — Va{b) are linear combinations of the coordinates 
of b with integer coefficients. The vertices of the simple polytope Pb are precisely 
those vectors which have all coordinates non-negative. The edges of Pb 
emanating from a vertex Va- are parallel to certain non-zero vectors with minimal 
support in the kernel of A. These vectors are called circuits in matroid theory. 
For any index i G {1, . . . ,n}\a let C{a,i) € { — 1,0,-1-1}" the associated basic 
circuit. This is the unique vector in the kernel of A whose support is a subset of 
aU{i} and whose z-th coordinate is +1. The following lemma is straightforward: 

Lemma 3.2 For the vertex of Pb indexed by a, the series ^ equals 

f{cone{Pb.v,),x) = x"-(^) • [](l - (7) 

In the formula above, only the monomial a;""^^-' depends on the specific right 
hand side vector b. The other factors depend only on the chamber which contains 
b. We use the following procedure for computing the generating function for the 
set of non-negative integer solutions to a unimodular system Ax = b. 

Algorithm 2. (Computing the BBKLP generating function) 

Input: Unimodular matrix A, a representative vector b for a chamber of A. 
Output: The generating function f{Pb, x) for the set of lattice points in Pb- 

1. Compute all the circuits of matrix A. This step is entirely independent of 
b and can be done a priori, before processing any particular chambers. 

2. List all subsets ct C {1, . . . , n} which index vertices of the polytope Pb- 

3. For each a in the previous step, compute the right hand side of (|^). 

4. Output the sum /(cone(Pb, Va), x) of these rational functions. 



We illustrate the output of this algorithm for the unimodular 3 x 6-matrix 

1 110 

-10 110 

-10-101 

and the right hand side vector b = 62, ^3) in the same chamber with (1, 3, —2). 
The polytope Pb is three-dimensional and has six vertices. Their index sets tr 
and corresponding generating functions /(cone(Pb, v^r), x) are listed in Table |l|. 
The lattice point enumerator f{Pb, x) is the sum of these six rational functions. 
The number of lattice points in Pb is found to be: 



/(n, (1,1, 1,1, 1,1)) 



libi + 2)(6i + l)(26i + 3&2 + 363 + 3). (8) 
6 
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index set 



rational function 



_ XiX4\ X3Xi \ f _ X4X6 \ 

\ X2 J \ X2X5 J \ Xs J 

( 1 _ ( 1 — £3_v^ ( I _ _^Y^ 

y xiXi J y XiX4,Xs j \^ X4Xe j 



{2, 4, 5} 
{1,4, 6} 
{2, 4, 6} 
{3, 4, 5} 
{3, 4, 6} 
{1,4, 5} 



X2^^X^~^^-^-^X^^^+^^+^-- 



2,1x42,6 

fcl m . 62 +63+62 { 1 2,5 \ ^ (-^ _ 2,12,4 \ "'^ ^ 2 — ^3 

2,42,6 y y 2-2 y \ 2,22,6 ^ 

6J . 62 ,T-_63+62 / 1 _ 2,22,6 A ^ 25 A ""^ 2^ _ 2,12426 ^ 

2,3 y \ 2,42,6 y \^ 2,3 



25 



2,1 24 



Table 1: Rational functions associated to the six supporting cones at the vertices 



This last evaluation can be done symbolically, for instance, using the com- 
mand simplify in Maple, but the symbolic simplification is too slow for larger 
examples. We compute the limit Xi ^ 1 in such rational functions by first 
specializing to a single variable t, in a manner to be described in Algorithm 3. 

We next discuss how to implement Step 2 of the algorithm, namely, how to 
efficiently list all vertices of Pb. The first possibility is to compute the prime 
decomposition of the monomial ideal M which was used in Section 2 to encode 
the chamber of b. Indeed, a subset a corresponds to a vertex of Pt if and only if 
{xi : i ^ a) is a minimal prime of M. The second possibility is to precompute 
the vector- valued linear functions Va-{b) for all column bases of A. Similarly 
to the computation of the circuits in step 1, this can be done a priori, before 
processing any particular chambers. For any particular chamber, we take the 
sum in step 4 only over those bases a which satisfy v„{b) > 0. In our practical 
implementation we opted for a third possibility, namely, to apply a depth-first 
search algorithm to the graph of basic feasible solutions of Ax = b,x > 0, 
where the edges are basis exchanges as in the simplex algorithm ||2^, Chapter 
8]. A considerable speed-up over our crude Maple implementation can still be 
obtained by using the reverse-search algorithm of Avis and Fukuda Q . 

The output of Algorithm 2 is a generating function /(Pb, x) that represents 
the vector partition function on a particular chamber. Now we face the prob- 
lem to evaluate, for any particular b G Z'', the hmit of /(Pb, a;) as x tends to 
(1,1, In the literature there are two approaches to this problem: the 
Barvinok-Brion method Algorithm 5.2] and the Dyer-Kannan method 
Both methods consider the rational series as a sum of exponential functions 
each of which converges for almost all choices of x. The first approach essen- 
tially takes the residue of the function and the second computes the value of 
the rational function at a point close to a; = (1, 1, . . . , 1) and carefully rounds 
the answer to the nearest integer. When we tried these two approaches ex- 
perimentally, we ran into memory problems and numerical instabilities. In our 
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experience, the following alternative method works rather well in practice: 
Algorithm 3. (Evaluating the BBKLP generating function at (1, 1,1,..., 1)) 

1. Eliminate rank{A) many variables by substitutions Xi = I where i runs 
over a column basis of A. All denominators 1 — a;*^^"^''^ remain nonzero. 

2. For each vertex Vu of Pb, replace each remaining variable xj by 1 — jt. 
This transforms (0) into a rational function in one variable t. We express 
the result in the form numerator/denominator, where the numerator and 
denominator are relatively prime polynomials in t with integer coefficients. 

3. Replace the sum of rational functions, one for each vertex of P^, by a 
single rational function p(t)/q{t). Here q{t) is the least common multiple 
of the denominators of the rational functions produced in step 2. 

4. Both p{t) and q{t) vanish at t — 0. Let be the largest common factor. 
We compute the limit of p{t)/q{t) as t ^ using L'Hopital's rule. For 
that we need the value at of the a-th derivatives of p and q. These can 
be found in Maple using the built-in feature of automatic differentiation. 
This allows us to retain the representation of p as sum of terms, one for 
each vertex of P^, and that of g as a product of binomials 1 — x'~'^'^'^\ 

5. Output (//5)(0). 

At the beginning of this section we had assumed that 6 lies in the relative 
interior of a maximal chamber. This assumption can be removed easily. It was 
made in order to uniquely identify the chamber and hence a representation of 
Pb as a simple polytope. If h happens to lie in a lower-dimensional chamber, and 
Pb is not simple, then we can use the combinatorial description of any adjacent 
maximal chamber in step 2 of Algorithm 2. This is consistent with the fact, 
implied by Theorem [l.l| , that the polynomials representing 4>A{b) on different 
chambers must agree on the intersection of the closures of these chambers. 



4 Contingency Tables 

In the remainder of this paper, we report on the implementation and perfor- 
mance of our methods for two important families of unimodular matrices A. 
We present both computational and mathematical results. We ran all our ex- 
periments in a computer with a single Pentium-Ill CPU with 700Mhz and 256 
MB RAM using the computer algebra packages Macaulay 2 and Maple. The 
generation of chambers was performed using Topcom and Puntos; see ^, |20| . 

Let r — (ri, . . . , r^) and c — (ci, . . . , c„) be compositions of a fixed integer 
A" > 1. Let Ere denote the set of all m x n non-negative integer matrices in 
which row i has sum and column j has sum Cj . Thus ^ Tij = N for 
any T G E^c- We are interested in the number fj^Y^rc of matrices in E^c- This 
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number equals 4>A{b) where A is the node-edge incidence matrix of the complete 
bipartite graph Kn^m and the vector b is the vector (r, c). Thus we are counting 
the lattice points in a transportation polytope. There is an extensive hterature 
on computing the function (r, c) t-^ ^Yir,c- See pO[ | and the references therein. 

We implemented the Grobner bases algorithm described in Section 2 in the 
computer algebra system Macaulay 2, which was developed by Grayson and 
Stillman fisll . Our Macaulay 2 program for computing the polynomial repre- 
senting (j)A on a single chamber is very short and simple. In Appendix 2 we list 
the entire program for one chamber in the 4 x 4-contingency table case. 

As mentioned in the introduction, 4 x 4-tables are an important benchmark. 
There are 3694 chambers modulo symmetry. On each chamber, the function 
(r, c) 4F^rc is a polynomial of degree nine in the eight variables t"!, 7*2, ?'3, ^4, 
Ci,C2,C3,C4. Mount computed (interpolation schemes for) all 3694 polyno- 
mials. He reported a 3 hour calculation for each chamber, adding up to a total 
of 6 weeks of distributed computing for preprocessing all chamber polynomials. 

Our experiments show that the Grobner basis computation is as least as fast 
as Mount's interpolation technique. We computed all 3694 chamber polyno- 
mials using the Macaulay 2 code listed in Appendix 2. The running time per 
chamber ranged from 7 seconds to 45 minutes. It took us 6 1/2 weeks sequential 
computing time to complete the task. Our Macaulay 2 code can easily be mod- 
ified to get the numerical value #Src for any given r, c G N^. Computing such 
numerical instances takes 20 seconds on the average for 4 x 4-tablcs. Similar 
computations for 4 x 5-tables have not yet been successful in Macaulay 2. 

We implemented the BBKLP method described in Section 3 for contingency 
tables in Maple. The generation of the BBKLP rational function (Algorithm 2) 
runs rather well for our purpose. It takes only a few seconds for 4 x 4-tables, 
as little as five minutes for 4 x 5-tables and up to two days for 5 x 5-tables. We 
wish to stress that our Maple code does not use optimal techniques for vertex 
enumeration of polytopes. For instance, using the Avis-Fukuda reverse search 
algorithm j|] instead of depth-first search would give a significant speed-up over 
our crude implementation. For example, the vertices of a 5 x 5 transportation 
polytope can be computed in a few seconds using the program Irs [Q. 

The second stage in the BBKLP method is Algorithm 3. This can be applied 
either for symbolic parameters Ti and Cj , in which case the output is a chamber 
polynomial, or for numerical values of and r^, in which case the output is 
the integer ^Tir^c- The second application of Algorithm 3 performs extremely 
well in Maple. The running time of a numerical evaluation (r, c) ^ 4i^'^r,c 
using Algorithm 3 is close to one minute for 4x4 tables, about ten minutes for 
4x5 tables, and about ten days for 5 x 5-tables. On the other hand, the first 
(symbolic) application of Algorithm 3 is only possible for smaller matrices, and 
is generally outperformed by the Grobner basis computation in Macaulay 2. 

Here are three test cases that show the power of the BBKLP technique, 
with numerical evaluation in Algorithm 3. The largest instance computed by 
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Mount is the number of 4 x 5-tables with margins [3046, 5173, 6116, 10928] 
and [182,778,3635,9558,11110]. It took him 20 minutes of parallel computing 
to find the value 23196436596128897574829611531938753. Our Maple program 
reproduces this number in only 10 minutes. 

Consider next the 4 x 5-tables whose margins are [338106, 574203,678876, 
1213008] and [20202, 142746,410755, 1007773, 1222717]. Their number equals 

316052820930116909459822049052149787748004963058022997262397. 

The computation took 35 minutes. The associated transportation polytope Pf, 
is 12-dimensional and has 976 vertices. 

Finally, we counted all 5 x 5-tables with margins [30201, 59791, 70017, 41731, 
58270] and [81016,68993,47000,43001,20000]. The associated 16-dimensional 
transportation polytope has 13150 vertices. This computation took 10 days and 
the answer is a 64 digit number. Algorithm 2 ran about 2 1/2 days. The size 
of its output exceeds the memory of our computer. Therefore we had to apply 
the Icm-computation in Algorithm 3 to incremental pieces of this output. 



Our Maple program for counting 4x4 and 4 x 5-tables is available at tfww . 



math.ucdavis.edu/~deloera/contingency.html. This webpage includes all 



relevant data for the two specific 4 x 5-tables discussed above. 

The subproblem of enumerating all chambers lead us to take a look at the 
structure of the chamber complex for the n x m contingency tables. This cham- 
ber complex is the cone over the chamber complex of the product of two simplices 

n m 

A„_ixA„_i = { (xi,...,a:„;2/i,...,y™) e M;i+™ : = 1}. 

i=i j=i 

The combinatorial structure of the polytope A„_i x A„i_i can be read off from 
the complete bipartite graph Kn,m- For instance, the full-dimensional simplices 
in A„„i X Am-i correspond to spanning trees of Kn^rm while the facets of 
A„_i X A„i_i are complete bipartite subgraphs of Kn^m obtained by removing 
a vertex. The (n + m — 3)-dimensional subsimplices correspond to a spanning 
tree minus an edge. We define a diagonal section of A„„i x Am-i to be any 
afhne hyperplane which is spanned by vertices of A„„i x Am-i but is not a 
facet hyperplane. The diagonal sections are in bijection with spanning forests 
of Kn^m which have exactly two components. Let ^n,m denote the subdivision 
of the polytope A„_i x Am-i defined by the diagonal sections. Equivalently, 
two points {x, y) and (x', y') in A„_i x A„i_i lie in the same open cell of ^n.m 
if and only if the lie on the same side of any hyperplane of the form 

a^ii H h Xj,. + H h yj^ = 0. 

We call iln,m the diagonal section complex oi A„ x A„i. 

Proposition 4.1 The chamber complex o/A„x A„j coincides with the diagonal 
section complex fln.rn- There exist virtual chambers whenever m + n > 7 . 
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Proof: For any polytope whatsoever, the diagonal section complex can be de- 
fined, and it always refines the chamber complex. The two complexes are equal 
for polygons, but they are usually not equal for higher dimensional polytopes. 
What we are claiming is that they are equal for products of two simplices. 

The key observation is this: the intersection of a diagonal section with any 
facet of the polytope A„ x Am equals the convex hull of all vertices of the facet 
which lie in that diagonal section. This follows from our graph-theoretical dic- 
tionary, since each facet corresponds to a complete bipartite subgraph Kn-i,m 
or Kn,m-i- From this it follows that each codimension one simplex spanned by 
vertices of A„ x A^ has the same intersection with the boundary of A„ x A™ as 
the corresponding diagonal section does. Therefore the chamber complex equals 
the diagonal section complex. The assertion about virtual chambers is proved 
by computer calculations for K2^5 and ^'3,4. | 



5 Kost ant's Partition Function 

Let A be the node- arc incidence matrix of the complete acyclic graph K^- 
The function (pK„ '■— 4>a is the Kostant partition function for the root sys- 
tem An-i- Explicitly, let ei, . . . ,e„ denote the standard basis of Z", and let 
-E'1,2, Ei,s, . . . , En^n~i denote the standard basis of zi^) . The matrix A repre- 
sents the map 

T : M.^''} R", E,^j ^e,- e-j for 1 < i < j < n. 
The image of r is the (71 — l)-dimensional cone 

Sn = {(wii ■ ■ • i^n) £ M" I • > for \ < i < n and ui+- ■ — O}. 

Kirillov [pT| page 57] posed the problem of finding the number of chambers 
for Kostant 's partition function. We give a partial solution to Kirillov's problem 
by determining the number of chambers for n < 7. See Table ^ below. 

We also computed all chamber polynomials representing (j)K„{b) for n < 6. 
This was done using our Macaulay 2 implementation (see Appendix 2) of the 
Grobner basis method in Section 2. For instance, for n = 6, there are 820 cham- 
ber polynomials, each of degree 10 in five variables. All of these polynomials 
are available, both in expanded form and as an on-line calculator, at our web 



site www .math .ucdavis . edu/~deloera/kostant .html 



As a small sample of our results we present all chambers and chamber poly- 
nomials for n = 4. These polynomials were first computed by mathematical 
physicists in pH |. Analogous computations for n > 5 had been infeasible in 
1984. In Appendix 1 we list all those chamber polynomials for n = 5 which 
can be factored over Q. Several authors |l] have studied factorization pat- 
terns of polynomials representing Kostant's partition function. A forthcoming 
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Figure 1: The chamber complex for the complete graph K^^ 



paper by Postnikov and Stanley contains the state of the art. Our data provide 
complementary information to their combinatorial results. 

The cone S4 spanned by the columns of the node-arc incidence matrix of 
is a three-dimensional triangular cone. The chamber complex is a subdivision 
of this cone into seven triangular cones. See Figure ^ for a 2-dimensional per- 
pendicular slice showing the chamber complex. The formulas below are given 
only in terms of 61, 62, ^3, in view of 64 = —61 — 62 — ^3- By the symmetry of the 
example it is enough to give the four polynomials for the indicated chambers in 
Figure |. The label of a chamber in the figure and its polynomial match. 

1. If min{fe3,-b2,bi+fe2} > then0K4(fe) = (6i-hfe2+3)(bi + 62 + 2)(6i-hb2 + l)/6. 

2. If min{6i,62,63} > then 0K4(&) = (61 -I- l)(&i 2)(bi 362 + 3)/6 

3. If min{6i,b2, 61 4-63,62 + 63,-63} > then (;/)K4(fe) = l + ^bi+2/3b3 + 
62+3/2 61 62 + 61" + 1/6 6i''' + 1/2 61^62 - 1/6 63^ -1/2 61 63" + 1/2 61 63 - 1/2 63^ 

4. If min{6i, 62+63, -61-63} > then 0^4(6) = (6i+2)(6i + l)(26i+362+3+363) 

Let T{Kn) be the chamber complex for (t>K„- This is a polyhedral decompo- 
sition of the cone S'„. We have the following result: 

Theorem 5.1 The complex T{Kn) has chambers with at least 2L"/^J facets. 
There exist virtual chambers for n > 5. The exact number of chambers for 
n < 7 is given by Table |. 

Proof: Let An-i ^ {ci ~ Cj : 1 < i < j < n}. There is a well-known bijection 
between cuts of the digraph Kn and hyperplanes spanned by subsets of An-i- 
For odd values of n there are "balanced" cuts for Kn- By a balanced cut we 
mean one where the hyperplane associated divides the set of roots — ej outside 
the hyperplane into equal size groups. In Figure |^ we show one such cut for 
K5 that leaves two roots in each side of the plane X3 = 0. To obtain such a 
balanced cut for general Kn, odd n, note that there is a middle node labeled 
[n/2j + 1 that has exactly as many entering arcs as leaving arcs. The cut 
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n 


Number of chambers 


Degree of (pK^ 


3 


2 


1 


4 


7 


3 


5 


48 


6 


6 


820 


10 


7 


44288 


15 



Table 2: Chambers for A. 




{1, 2, 3, . . . , [n/2j , [n/2j + 2, . . . , n} and { [n/2j + 1} is balanced. The vectors 
in A„_i that lie on the plane a;[ra/2j+i = form the configuration An-2- 

The intersection of the chamber complex of A„_i with a balanced hyperplane 
H induces exactly the chamber complex of A„-2- Indeed, the only way to create 
new cells for H n An-i (not already in An-2) is if simplices with vertices on 
opposite halfspaces of H cut out new vertices in H . But pairs of vectors on 
opposite sides of H are always collinear with a root Ci — ej lying on H . The 
collinearities can be read off from cycles of length three in the graph Kn that 
touch the vertex [n/2j + 1. The existence of triples of collinear vectors, the 
center one inside H, has another effect: a chamber 7 of i?, one of whose vertices 
is part of a coUinearity, extends to both halfspaces of H. This is because the 
(n — 2)-simplices inside H that make up that chamber can be turned into (n — 1)- 
dimensional simplices by coning them with the two extremes of the collinearity 
that do not belong to the hyperplane H . Note that the completion happens in 
both halfspaces of H but the result of intersecting these simplices has in common 
the open cell 7 that connects both sides. This might not be the final chamber 
that extends 7, as other vectors in An-i not lying on H could be used to build 
and intersect more (n — 1) simplices, but the the result will be contained in this 
initial convex cell that touches both halfspaces of H. The number of facets will 
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be then at least twice the number of facets of 7. The doubling on the number of 
facets occurs for odd values of n but for even values at worse remains the same. 
Thus, recursively we can build a chamber with exponentially many facets. The 
rest of the statement follows from computer calculations based on the duality 
between chambers and triangulations as explained in Section 2.| 
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6 Appendix: Kostant partition function for A4 

Here we consider (px^ for n = 5. This is Kostant's partition function for the 
root system A4. The chamber complex can be visualized as a subdivision of a 
tetrahedron. This polyhedral complex has 19 vertices, 77 edges, 107 triangles 
and 48 three-dimensional chambers. Only two of these 48 chambers are not 
tetrahedra: they are bipyramids. Thirty of the 48 chamber polynomials are 
irreducible over Q. Wc explicitly list the other 18 chamber polynomials, namely 
those that factor, together with defining inequalities for their chambers. 

1. If min {&i, ?)2, 63, ^4} > then 

(61 + 3) (61 + 2) (61 + 1) (61^ + 56162 + 96i + 20 + 106i + 3O62) 
(62 + 3 + 61 + 363) 

2. If min {64, 61 + 62 + 63, -63, -62} > then 

glo (62 + 5 + 61 + 63) (62 + 4 + 61 + 63) (62 + 3 + 61 + 63) (62 + 2 + «1 + 63) 
(62 + 1 + 61 + 63) (62 + 3 + 61 - 263) 

3. If min {—63, —62 — 64, —61 — 64, 61 + 6a + 63 + 64} > then 

gig (64 + 3 + 61 + 63 + 62) (64 + 2 + 61 + 63 + 62) (64 + 1 + 61 + 63 + 62) 

(60 + 566i + 662 - 1463 - 5464 + 9636461 + 66263 - 3626461 - %ibi + 

36i6| + 3626I - 6636I + 276162 - 96163 - 96?64 + 66i64 + 96?62 - 661 6i + 36i64 + 

66362 + 246364 - 456164 - 6626364 - 6i + 66? - 96i - 1563 - 26f + 36i + 96l) 

4. If min {-63, -62, -64, 61 + 62 + 63 + 6^ } > then 

3^0 (62 + 3 + 61 - 263) (64 + 3 + 61 + 63 + 62) (64 + 2 + 6j + 63 + 62) 

(64 + 1 + 61 + 63 + 6g)(6i + 26263 + 962 + 26162 - 36264 + 963 + 6? + 20 + 96i + 

26163 + 6i - 36364 - 2I64 - 36164 + 66l) 

5. If min {63, 64, —62, 61 + 62} > then 

3^0 (62 + 61 + 3) (62 + 2 + 61) (62 + 1 + 61) (62 + 5 + 61) (62 + 4 + 61) 
(62 + 3 + 61 + 363) 

6. If min {61, 63, 62 + 64, —61 — 63 — 64} > then 

^ (61 + 3) (61 + 2) (61 + 1) 

(60 + 566i + IIO62 + 7O63 + 5O64 - 3O636461 + 9O6263 + 306^63 - I562646J 
-66^63 - 156361 - 306164^ - 3O626I - 3O636I + 576162 + 2I61 63 - 156? b4 + 36? 62 + 
156i6i - 3063^64 + 3O6264 - 156164 + 15626361 - 106i - 206t + 661^ + 606i - 26f + 
106i - 306|) 
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7. If min {61, 63 + 64, 62 + 64, — &i — 64} > then 

+ (61+2) (61 + 1) 
{3b1b3 - 661^64 + 9&? + 36?62 + 516i + 576162 + 156i6i - 96164 
-156i6l + 276163 - 15626461 + 15626361 + 3O6264 + 60 + 6O63 + 606i + 4O64 + 
9O6263 - 306| + IIO62 - 106t + 3062^63 + 106i - 3O626I) 

8. If min {—62 — 64, —62 — 63, —61 — 63 — 6^ , 61 + 62 + 63 + 64} > then 

(64 + 3 + 61 + 63 + 62) (64 + 2 + 61 + 63 + 62) (64 + 1 + 61 + 63 + 62) 
(60 + 5I61 + II62 - 963 - 5964 + 3636461 - 66^63 - 9626461 - 36?63 
-66361 +96161-96364^+276162 -96163 -361^64 + 96i64 + 66^62 -36161 + 66364 + 
24636^ - 396164 + 6626361 - 363-'' + 26l + 96? - 126i - 1863^ + 6i + 126^) 

9. If min {—62, —63 — 64, 61 + 62 + 63 + 6^ , —61 — 62 — 64} > then 

(64 + 3 + 61 + 63 + 62) (64 + 2 + 61 + 63 + 62) (64 + 1 + 61 + 63 + 62) 
(62+4-64-63+61) {2bl + 1362 + 46162 - 6264 - 6263 + 15 + 136i 
-6163 + 26? - 763 - 764 - 6164 + 26| + 46364 + 26i) 

10. If min {64, 61 +62 + 63,-62 - 63, -bi - 63} > then 

- 3I0 (^2 - 3 + 63 - 26i) (62 + 5 + 61 + 63) (62 + 4 + 61 + 63) 
(62 + 3 + 61 + 63) (62 + 2 + 61 + 63) (62 + 1 + 61 + 63) 

11. If min {61, 62 + 63 + 64, —61 — 64, —hi — 63} > then 

^ (61 + 3) (61 + 2) (61 + 1) (63 - 264 + 62+3) (106i + 3O62 + 156162 + 2O6263 + 
2O6264 + 20 + 156164 + 3O63 + 3O64 + 156163 + 66? + 1064^ + 1061 + 246i + 2O6364) 

12. If min {61, 64, 62 + 63, -61 - 63} > then 

3^0 (&i + 3) (61 + 2) (61 + 1) (61 + 4 + 262 + 263) 

(56i + 56162 + IO6263 + 2O62 + 136i + 56163 + 2O63 + 26? + 56i + 15) 

13. If min {-62 - 63, -64,61 + 62 + 63 + 64,-61 - 63} > then 

-3I0 (&2 - 3 + 63 - 26i) (64 + 3 + 61 + 63 + 62) (64 + 2 + 61 + 63 + 62) 
(64 + 1 + 61 + 63 + 62) (6i + 26263 + 962 + 26162 - 36264 + 963 
+6? + 20 + 96i + 26163 + 6i - 36364 - 2I64 - 36164 + 66l) 

14. If min {-62, 63 + 64, 61+62, -61 - 62 - 64} > then 
3^0 (62 + 1 + 61) (62 + 5 + 61) (62 + 4 + 61) (62 + 61+3) 
(62 + 2 + 6i)(262 + 26i + 364 + 3 + 363) 

15. If min {—61 — 64, —62 — 63 — 64, 61 + 62 + 63 + 64, —61 — 63} > then 
3ljy (63 - 264 + 62 + 3) (64 + 1 + 61 + 63 + 62) (64 + 3 + 6j + 63 + 62) 
(64 + 2 + 61 + 63 + 62) (6? + 26264 + 26263 

-662 - 36162 + 20 + 6l + 26364 + 6i - 664 - 36164 - 66s - 36163 + 246i + 66?) 

16. If min {61, 62, 63 + 64, —61 — 62 — 64} > then 

^ (61 + 3) (61 + 2) (61 + 1) (6? + 56162 + 96i + 20 + 106i + 3O62) 
(262 + 26i + 364 + 3 + 363) 

17. If min {61, 62 + 63 + 64, —61 — 63 — 64, —61 — 62 — 64} > then 

^ (61 + 3) (61 + 2) (61 + 1) (-64 + 3 - 63 + 262) (1061+3062 + 156162+206263 + 
2O6264 + 20 + 156164 + 3O63 + 3O64 + 156163 + 66? + 1061 + 106i + 246i + 2O63 6^ ) 

18. If min {—62 — 63 — 64, —61 — 63 — 64, bi + 62 + 63 + 64, —61 — 62 — 64} > then 

(-64 + 3 - 63 + 262) (64 + 2 + 61 + 63 + 62) (64 + 1 + 61 + 63 + 62) 
(64 + 3 + 61 + 63 + 62) (6? + 26264 + 26263 - 662 - 36162 +20 + 61 
+26364 + 6i - 664 - 36164 - 663 - 36163 + 246i + 66?) 
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7 Appendix: Macaulay 2 program 



In this appendix we present our implementation of the Grobner basis algorithm 
from Section 2. For an introduction to the computer algebra system Macaulay 
2 see and fl^ . Our program starts by defining the unimodular 8 x 16-matrix 



A of rank 7 which 


represents the counting 
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We next input a weight vector W of length 16, to be interpreted as a 4 x 4-table: 
W = {1000,1,1,1, 1,1,1,1, 1,1,1,1, 1,1,1,1} 

The following five command lines compute the monomial ideal M = m^, {Ja)i 
here called nonf aces, which represents the chamber we are interested in: 

n=#W; d=n-#A+l;R= qq[x_l. .x_n, Weights => W] 
Binomial = (b,R) -> ( top := 1_R; bottom := 1_R; 

scan(#b, i -> if b_i > then top = top * R_i~(b_i) 

else bottom = bottom * R_i~(-b_i)); top - bottom); 
nonfaces = ideal leadTerm ideal applyC A, a -> Binomial(a,R) ) ; 

We compute the presentation ideal M + La of the cohomology ring H*{Xy^] k). 
It is denoted 1. 

S = QQ[x_l..x_n, rl,r2,r3,r4,cl,c2,c3,c4] ; 
f = map(S,R, toList (x_l . .x_n) ) ; 

Linform = (b,R) -> (s := 0_R; scan(#b,i -> s = s + b_i*R_i) ; s) ; 
I = f (nonfaces) + 

ideal apply (entries transpose syz matrix A, a -> Linf orm(a, S) ) ; 

The next four lines compute a representation of the Todd class modulo 1 . 

todd = (x) -> (l+l/2*x+l/12*x-2-l/720*x-4+l/30240*x-6-l/1209600*x-8) ; 
trunc = (d,f) -> sum select(terms f, t -> sum degree t < d+1) ; 
toddclass := 1_S; 

scan(l..n, i -> toddclass = trunc(d, toddclass * todd(x_i)) 7. 1); 

All subsequent computations take place in the quotient ring T = S/l = 
H*{Xw- k). We compute all successive powers of a general divisor ^ UiXi. 
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T = S/I; 

g = map(T,T,joiii(toList(ii:l) , {rl ,r2,r3,r4, cl,c2,c3,c4})) ; 
u= (0, rl-c2-c3-c4,c2,c3,c4,r2,0,0,0,r3,0,0,0,r4,0,0,0) ; 

divp = 1 ; 

divpowers = apply (l..d, i -> 

(divp = sum toList apply (l..n, i -> u_i * x_i * divp))); 

In the final four lines of code, the graded components of the Todd class are 
multiplied with the complementary powers of the divisor '^UiXi. The products 
are added up (in T) and the sum is normalized so that its constant term is 1: 

component = (d,f) -> sum select(terms f, t -> d == sum degree t) ; 
erhart = sum toList apply (0 .. d-1 , 

i -> (divpowers_i * (l/(i+l)!) * component (d-i-l,toddclass)) ) ; 
toString (g(erhart)/g(component(d,toddclass)) + 1) 

The final output is a polynomial of degree 9 in the variables rl, r2, r3, r4, cl, 
c2, c3, c4. This particular chamber polynomial has 1967 terms. The running 
time of this entire piece of code is about 25 minutes. 

Users of Macaulay 2 will find it easy to modify our code so that it works for 
any unimodular matrix A and any right hand side b = Aw. Besides redefining 
the variables A and w, one only needs to change those command lines which 
involve the variables rl, r2, r3, r4, cl, c2, c3, c4 particular to 4 x 4-tables. 
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